#pragma once

#include "vector.h"
#include "smatrix.h"
#include <cassert>
#include <iterator>

namespace g2 {

template <typename T>
class basic_line;

typedef basic_line<double> line;
typedef basic_line<float> linef;
typedef basic_line<long double> linel;

/// represents the line on plane: a*x + b*y + c = 0
template <typename T>
class basic_line {
public:
    T a, b, c;
public:
	basic_line()
		{ }
    basic_line(T _a, T _b, T _c) 
		: a(_a), b(_b), c(_c) { }
	template <class U> 
	basic_line(const basic_line<U> & v) 
		: a(T(v.a)), b(T(v.b)), c(T(v.c)) { }
public:
	basic_vector<T> projection(basic_vector<T> const & p) const;
	/// builds line that best fits given range of points in least squares sense
	template <typename I>
	void best_fit(I begin, I end);
};

template <typename T>
inline bool operator == (const basic_line<T> & l, const basic_line<T> & r)
{
	if (l.c == 0)
		return r.c == 0 && l.a == r.a && l.b == r.b;
	if (r.c == 0)
		return false;
	return
		r.c * l.a == l.c * r.a && r.c * l.b == l.c * r.b;
}

template <typename T>
inline bool operator != (const basic_line<T> & l, const basic_line<T> & r)
	{ return !(l == r); }

template <typename T>
inline basic_vector<T> basic_line<T>::projection(basic_vector<T> const & p) const
{
	basic_vector<T> const n(a, b);
	return p - (c + scalar_product(p, n)) / sqr(n) * n;
}

template <typename T>
template <typename I>
void basic_line<T>::best_fit(I begin, I end) 
{
	assert(std::distance(begin, end) >= 2);
	basic_vector<T> sum1(0, 0);
	basic_smatrix<T> sum2(0, 0, 0);
	for (; begin != end; ++begin) {
		basic_vector<T> const i(*begin);
		sum1 += i;
		sum2 += tsquare(i);
	}
	if (sum2) {
		basic_vector<T> ab = sum2.inverse(sum1);
		a = ab.x;
		b = ab.y;
		c = -1;
		return;
	}
	c = 0;
	if (sum2.yy != 0) {
		a = 1;
		b = - sum2.xy / sum2.yy;
		return;
	}
	a = - sum2.xy / sum2.xx;
	b = 1;
}

} // g2
